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Abstract 

We introduce a new transfer matrix method for calculating the thermody- 
namic properties of random-tiling models of quasicrystals in any number of 
dimensions, and describe how it may be used to calculate the phason elastic 
properties of these models, which are related to experimental measurables 
such as phason Debye- Waller factors, and diffuse scattering wings near Bragg 
peaks. We apply our method to the canonical-cell model of the icosahedral 
phase, making use of results from a previously-presented calculation in which 
the possible structures for this model under specific periodic boundary condi- 
tions were cataloged using a computational technique. We give results for the 
configurational entropy density and the two fundamental elastic constants 
for a range of system sizes. The method is general enough allow a similar 
calculation to be performed for any other random tiling model. 
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I. INTRODUCTION 



Quasicrystals are defined as structures which possess translational order to the extent 
that their Fourier transforms exhibit 5-function Bragg peaks, but which have symmetries 
that are forbidden in a periodic Bravais lattice. A number of alloys have been found ex- 
perimentally which appear to be true quasicrystals in this sense, such as i(AlCuFe) [I[ and 
i(AlPdMn) PJ^j, which both exhibit resolution-limited Bragg peaks. 

Two competing physical scenarios have been advanced to explain the origins of quasiperi- 
odic ordering. One of them — the 'ideal tiling' scenario — postulates that the atomic structure 
of a well-annealed quasicrystalline sample is perfectly quasiperiodic, like the two-dimensional 
Penrose tiling Q or its generalizations. In this scenario it is hypothesized that the micro- 
scopic Hamiltonian constrains the local arrangements of atoms in such a way as to implement 
something akin to Penrose's local 'matching rules', which force long-range quasiperiodicity. 

The alternative to this ideal tiling approach is the 'random tiling' scenario, which also 
makes use of local clusters of atoms, with packing rules similar qualitatively to Penrose's 
rules but insufficiently strong to force a unique behavior at long distances. This gives rise 
to an ensemble of different ways of packing atoms into the space occupied by the sample. 
Normally we assume these different configurations to be nearly degenerate in energy. In 
the actual models studied, space is assumed to be filled by a finite set of local patterns 
which we represent by a set of tiles. Such models have a contribution to their entropy 
arising from the many ways in which the tiles may be packed. This entropy can reduce 
the free energy with respect to crystalline structures (which are expected to be stable at 
zero temperature), offering an alternative explanation of how the quasicrystal state might 
be stabilized |5|||. To date, random-tiling models have been studied in two dimensions 
with 8-fold ]7|||, 10-fold [||,[|[10| and 12- fold symmetries fill , and in three dimensions with 
icosahedral symmetry |12]-[14|. 

In principle it should be possible to construct random-tiling models so that there is a 
one-to-one mapping of tiling configurations onto real atomic configurations, via deterministic 
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rules that specify how the tiles are to be decorated with atoms. Then the configurational 
entropy (and other thermodynamic quantities) of the atomistic model will be identical to 
that of the tiling model. The decoration of random tilings to produce atomic models of the 



icosahedral phases is discussed in Refs. [jT5-IB 



This paper is principally concerned with random tiling models of quasicrystals. The 
qualitative predictions of these models and the ideal tiling models are very similar, to the 
extent that they are known. In particular both models predict the sharp Bragg peaks 
observed in experiments. Experiments aimed at resolving the differences between the two 
for real quasicrystals have still not conclusively settled the question one way or the other. 
However, there are some arguments we can present in favor of the random-tilings. One 
important point is that we have some understanding of how interatomic potentials might 
lead to a random-tiling structure. The simplest model for such a process is the 2D 'binary 
tiling', in which atoms of two different sizes aggregate in patterns which can be set in one- 
to-one correspondence with the configurations of a tiling of 'fat' and 'skinny' rhombi j9]. 
More recently, it has been shown in simulations that several 3D models with reasonable 
pair potentials for simple atoms (Al, Mg, Li, Zn) will freeze into structures which can be 
represented as configurations of a random tiling model flT9| , |58| ]. On the other hand, there 
is no understanding of how the thermodynamics of a system governed by a microscopic 
Hamiltonian might give rise to the matching rules necessary for the ideal tiling scenario. In 
fact, the most promising recent advance in this direction has been achieved at the expense of 
blurring the distinction between random tilings and ideal quasicrystal models: one assumes 
an ensemble in which any of a large number of tilings is permitted, but adopts a Hamiltonian 
such that the ground state, and the thermodynamic state if the temperature is not too large, 
is essentially an ideal quasicrystal [EUJ. 
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A. Phason elasticity 



The long-wavelength behavior of random tiling models is described by a sort of Landau 
theory, which has two kinds of parameters: (i) the entropy density, and (ii) the 'phason 
strain' (see Sec. II). Most of the proposed experimental tests of the random-tiling scenario 
hinge around the expected variation of the entropy with this phason strain, which produces 
a kind of elastic term in the free energy. We describe this in more detail in the next section 
where we define the 'phason elasticity tensor', which measures the strength of this effect. 
It is predicted that the smallest eigenvalue of the phason elasticity tensor should decrease 
with decreasing temperature, producing instabilities of the quasicrystal phase with respect 
to other structures when the temperature becomes low enough |]2ll-p3|. 

Also, by contrast with the matching rule models, it is predicted that the phason elasticity 
for a random-tiling should have a gradient-squared form (see Section II A). This leads to 
diffuse 'wings' around peaks in the diffraction pattern, and to phason contributions to the 
Debye- Waller reduction of the Bragg intensities, both of which can in principle be measured 
experimentally to determine values for the elastic constants. 

B. Random canonical-cell tilings 

In this paper, we build upon results from our previously published work [J24|J , henceforth 
referred to as Paper I, to calculate a variety of properties of a random tiling ensemble 
based on the 'canonical cell' model of a quasicrystal. This model has been described in 
detail by Henley |25j . Briefly, it consists of vertices or 'nodes' which represent icosahedrally- 
symmetric clusters of atoms. Nearest-neighbor nodes are joined by 'linkages' of two types. 
Type b linkages run parallel to the axes of two-fold symmetry of the reference icosahedron 
and are all of the same length which we denote b; type c linkages run parallel to the axes of 
three-fold symmetry, and all have length c = by/3/2. The linkages can form three kinds of 
polygon which in turn form the faces of four kinds of tiles or 'canonical cells' A, B, C, and 
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D into which the entire space is divided. (These cells were depicted in Fig. 1 of Paper I and 
also Fig. 1 of Ref. [0.) 

In Paper I we gave an algorithm for generating a so-called 'stacking graph' for a random 
tiling, which we applied to the particular case of the canonical-cell tiling. This graph contains 
information about all possible structures that can be built out of a certain set of tiles with 
specified periodic boundary conditions. In this paper, we use this stacking graph as the basis 
for a transfer matrix calculation of the thermodynamic properties of the ensemble of random 
tilings of canonical cells, including the random tiling entropy and the phason elasticity J26J . 
These quantities have been found previously for another icosahedral random tiling, that of 
the two Ammann rhombohedra ||14|| . The canonical-cell model, however, is different in at 



least two respects: 

(i) It lends itself to the construction of realistic atomic models for icosahedral quasicrystals 
such as i(AlCuLi) 0, i(AMnSi) |TJ, i(TiCrSi) 0, and potentially i(AlPdMn). 



It is possible to construct decorations of the canonical cells which incorporate our 
understanding of the basic atomic motifs found in these alloys and which contain no 
points of unrealistically close or loose packing. Conversely it has not proved possible 
to relate the rhombohedron tiling to any good structure model. 

(ii) The canonical cell model is, unfortunately, much less tractable technically than the 
tiling of Ammann rhombohedra. In particular, Monte Carlo simulation has not been 
practicable because there is no known move that rearranges tiles locally and satisfies 



ergodicity. A possible chain-like or cluster move has been suggested by Oxborrow |]T5 



but this move is still not completely understood. (It has been understood in a two- 



dimensional toy model, the square-triangle random tiling [|TTJ, which shares some fea- 
tures with the canonical-cell tiling.) Since we suspect that this kind of intractability 
is more generic than the behavior of rhombus and rhombohedron tilings |Tl|j27| , |28[] , it 
behooves us to develop appropriate methods for studying less agreeable tilings such as 
the canonical-cell tiling, even if some other tiling ultimately proves more relevant to 



real quasicrystalline phases. 



C. Outline of the paper 

The paper is organized as follows. In Section II we develop the elastic theory of random 
tiling ensembles, including the special cases that apply to the ensembles we will be studying 
with our transfer matrix method. In Section III we define our transfer matrix and explain 
how it is used to calculate thermodynamic properties of the ensemble as a function of phason 
strain. In Section IV we describe the calculations we have performed using the canonical-cell 
tiling, and give our results for the configurational entropy and phason elasticity of the tiling 
for a variety of system sizes. In Section V we present our conclusions. 



II. PHASON ELASTICITY IN RANDOM TILINGS 

In studying icosahedral tilings, it is convenient to use a basis of unit vectors ejj, pointing 
from the center to vertices of the reference icosahedron, i.e., along the fivefold symmetry 
axes. It turns out that every type b linkage can be written as the sum of four of these basis 
vectors and every type c linkage can be written as the sum of three of them. Thus, having 
arbitrarily assigned r = at one node, every other node in a canonical-cell tiling may be 
represented in the form 
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r=^n Q e|, (1) 

a=l 

where n a are integers. The {ejj,} are independent over integers so this representation is 
unique. 

We can also define a 'phason' coordinate for each node 

r^Evi. (2) 

a=l 

The vector r 1 - lives in a space called 'phason' space or 'perp' space. As discussed else- 
where [21], with a proper choice of basis vectors {e^;} the perp-space coordinate is well- 



defined. We can view Equations (|l]) and @ together as defining a 'lifting' of the nodes to 
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points on a six-dimensional hypercubic lattice; thus an arbitrary tiling configuration may be 



viewed as the projection of a three-dimensional surface embedded in six-space |[29|| . Following 
the conventions of Jaric [BOH, we write 



e"i = 77(r,0,l), e ll 2 = 77(r,0,-l), e 1 ^ = 77(1, r, 0), 
ej[ = 77(0,1,7-), ef = 77(0,-1, r), e| = 77(1, -r, 0), 



(3) 



and 



e\ = 77(1, 0, -r), eS = 77(1, 0, r), = ^(-r, 1, 0), 

e ± 4 = 77(0, -r, 1), e-S, = 77(0,7-,!), e\ = ^(-r, -1, 0), 



(4) 



where r = |(1 + a/5) and 77 = (1 + r 2 ) ■ L / 2 |3l|. With these definitions, the length of the 



6-type linkage is 6 = 2y 1 + 2/y/E. We also define a coarse-grained phason coordinate 



h(r) = (r^), 



(5) 



where (. . .) r is an average over a neighborhood of r. The phason strain tensor is then the 
gradient of h in real space: 



Efr) = V P h(r). 



(6) 



A. Form of the elasticity tensor 

Consider an ensemble of tilings, possibly weighted with Boltzmann probabilities derived 
from some Hamiltonian 7i. We hypothesize that the Hamiltonian is sufficiently small, or the 
temperature sufficiently high, that difference in energy between configurations is typically 
much less than kT, so that the free energy of the ensemble will be dominated by the entropy. 
We further make the assumption that the relative weighting of long-wavelength phason 
fluctuations (after coarse-graining) will be determined by a free energy of the form [p!^ , |33|l 

T=\ rf 3 r/(E(r)) (7) 
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where 

/(E) = / + ^E-K-E + O(E 3 ). (8) 

K is the elasticity tensor (really a 4-tensor) and V is the system volume. In the case where 
the Hamiltonian TC is zero for all states in the ensemble, fo is equal to minus the entropy 



per unit volume Sy ||34j| . It is important to be aware that this form for the free energy is 
just a hypothesis. We cannot prove that /(E) must have an analytic minimum at E = 0, 
though it can be shown that, if the linkages in our tiling are uniformly distributed over 
all orientations equivalent by icosahedral symmetry, then r 1 - is approximately constant. In 
other words, the hypersurface formed by the tiling in six-dimensional space is approximately 
flat, having only small phason fluctuations at non-zero wavevectors. Thus the phason strain 
can be thought of as parametrizing the deviation from icosahedral symmetry, making it 
natural for / to have a minimum at E = if the free energy is dominated by the entropy 
of the ensemble. It would in theory be possible for / to have a non-analytic minimum at 
E = 0. However, in all models previously studied the analytic form (|J) as been borne out 
both by analytic calculations and by simulations p^JT3|j35[1 . Furthermore, a random tiling 
has long-range order (i.e., Bragg peaks exist) if and only if h(r) has finite variance; this can 
be shown to follow from the gradient-squared form (§]). 

It is common to deal with (H) by assuming that the net phason strain is zero and then 
re-expressing in terms of Fourier components h(q) (this generalizes trivially to the case 
of any fixed background phason strain — see Section II B below). The result is a sum over 
wavevectors q of terms each of which is a quadratic form in the three components of h(q). 
The coefficients are linear in {K{\ and quadratic in the components of q. This form is 
simpler than ([8]), since we write a quadratic form in three instead of nine quantities (and, 
what is more, the different q vectors decouple), and is appropriate for analyzing fluctuations 
from Monte Carlo simulations. However, the ensemble we will be working with in the present 
study permits undetermined phason strains with respect to only the (real space) z direction, 
and we measure only the elasticities associated with that direction. In this case the Fourier 
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analysis is not useful, and we must struggle with the full form of the elastic free energy. 



The elastic free energy for icosahedral symmetry can be written p0| , |36| , p7 | 



/(E) = \k x £ /•,/•,, + \k 2 {{E u + E 22 + E 33 ) 2 - | J] EaEv + 

(tE 12 + t~ 1 E 21 ) 2 + cyclic permutations} + 0(E 3 ). (9) 

If the icosahedral phase is to be stable against decomposition into other, possibly periodic 
phases, the quadratic stationary point in the free energy must be a minimum, rather than 
a saddle point or maximum. This requires that K\ be positive definite, and K 2 may have 
either sign but must lie within the range 

- ~ < K*/K x < ~ (10) 

The values of K\ and K 2 are one of the main points of contact between quasicrystal ex- 
periments and random-tiling theory. Knowing merely the ratio K 2 /K\ between the two is 
enough to determine the form of the diffuse scattering from the random tiling. 



B. Background phason strain 

Consider now the case of a tiling with periodic boundary conditions in all three real space 
directions. The periodicity of the structure is represented by a set of vectors {a,}, which are 
the displacement in real space from a node to the corresponding node in another periodic 
repetition of the structure. Since the vectors a^ connect nodes, they can be written in the 
form (]]]), each with a corresponding perp displacement a^j. The boundary conditions then 
constrain the system to have a global average phason strain B, defined by the linear system 
of equations 

B.a i = a J - i . (11) 

It will still be possible to consider fluctuations in the phason strain at finite wavevectors 
around this state, but when we consider the free energy associated with the phason elastic- 
ity, we must allow for the constraint on the long-wavelength phason strain by introducing 
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Lagrange multipliers. We will not follow through the complete analysis here, but merely 
quote the important results, which are 

(i) The elastic free energy density / is shifted by a constant term /q = |B ■ K ■ B. 

(ii) The quadratic term becomes ~(E — B) • K • (E — B). 

(iii) Most significantly, cubic and higher anharmonic terms now give contributions to the 
new quadratic term. Thus the elasticity tensor in the presence of the background 
phason strain loses its icosahedral symmetry and retains only the symmetry of the 
Bravais lattice defined by the periodic boundary conditions. The deviations, however, 
will be at most of 0(B) (from the cubic terms). Such terms in the elasticity were first 
measured by Oxborrow [|TT| . 

C. Elasticity theory for a tower 

In most of this paper, we consider a 'tower' in which periodic boundary conditions are 
applied in only the x and y directions; in the z direction the tower can be arbitrarily large. 
The towers we will be considering have a square base and, of the icosahedral symmetry 
operations, only x, y, and z reflections will be symmetries of the ensemble of possible 
configurations of cells in the tower; in particular, notice that the x and y axes will be 
inequivalent, even though the dimensions of the base in the x and y directions are equal, 
since the icosahedral point group has no 4-fold axis of symmetry. 

For such a tower, only gradients with respect to the real-space z coordinate remain free. 
Gradients with respect to x and y are fixed by the periodic boundary conditions, giving rise 
to a background phason strain of the type discussed above. We define B for a tower just as 
before except that the components B a3 are zero. Then 

B u = B 22 = B (12) 

with 
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£? = r -(2n+3) ; (13) 

where n is the order of the approximant, with n — 1 for the | structures, n = 2 for the 
|, and so forth. All other components of B are zero. For the phason strain tensor E, the 
boundary conditions give us just three free components E a % out of the nine. These three 
components transform like a vector, which we will call A, with 

A a = dh a /dz = E a3 . (14) 

We can write the free energy density for the tower in terms of this vector as 

/(A) = fo + l {CiA\ + C 2 A\ + C 3 (A 3 - 4 0) ) 2 } . (15) 

There are no cross-terms, because of the reflection symmetry, but the linear term shown 
involving the constant A^ is possible because the z reflection symmetry reverses both h 3 
and z. 

Using Equation ([J), we can show that the constants in this formula are given by 

C 1 = K 1 + (r 2 - ^)K 2 
C 2 = K X + (r- 2 - ^)K 2 

C 3 = Kt - l -K 2 (16) 

and 

4 0) = -1BK 2 j{K x - l -K 2 ). (17) 

We will calculate values for C%, C 2 , C3, and A3 using a transfer matrix method. Since 
these quantities are all functions of two elastic constants Ki,K 2 , there is some redundancy 
in such a calculation, which allows us a check on the accuracy of our methods. Inaccuracies 
can be ascribed to either (i) finite-size effects or (ii) finite- B effects (through the cubic terms 
in the elasticity theory, which we omitted). 
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D. Frequencies of objects in tiling 



The relative frequency with which various types of tiles appear in a random tiling is 
directly related to the phason strain. For example, in the case of certain extreme phason 
strains, the tiling can be composed entirely of one type of tile. In some random tilings, 
the number densities of the different tile types is fixed exactly by the phason strain. The 
two-dimensional Penrose tiling of rhombi is one such example. In the canonical cell tiling 
this is not the case, but there are still things we can say about the relationship of the number 
density of tiles to the phason strain. These topics were discussed in detail in Ref. P5|. Here, 
we will just summarize the situation. 

In addition to the phason strain, we can define many other macroscopic variables as the 
density (per unit volume) of occurrences of selected local patterns in the tiling. There are 
two reasons why we are interested in such patterns: 

(i) When we come to decorate our tilings with real atoms to create possible atomic struc- 
tures for icosahedral alloys, the simplest choice we can make is to use a single decoration 
for each type of canonical cell. In some cases a more satisfactory structure can be pro- 
duced using a "context-dependent" rule, in which the decoration of a tile depends on 
the those of the tiles surrounding it. In either case the decoration depends on the local 
patterns in the tiling and the density and stoichiometry of the resulting structure are 
simple functions of the densities of those patterns. 

(ii) As discussed below in Section III D, iteration of the transfer-matrix always generates 
an ensemble which maximizes the entropy per layer, whereas we actually want the 
ensemble which maximizes the entropy per unit volume. We will be introducing chem- 
ical potentials which couple to the densities of local patterns in the tiling. By varying 
these chemical potentials, we can vary the relative frequency in the tiling of the various 



patterns and so generate an ensemble which is closer to the desired one |38| . 
The simplest and most important of the density variables are the densities of the tiles 
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A, B, C, and D. In fact, these densities can be parameterized by just one independent 
parameter (, which is defined as the volume ratio V(D)/[V(A) + V(D)], where V(X) is the 
total volume occupied by tiles of type X. (This follows because V(B) and V(C) are unique 
functions of the phason strain — see Eqs. (3.1) and (3.6) of Ref. ||25| .) This is important not 
only because ( determines the stoichiometry and density of simple decoration models, but 
also because it completely determines the density of nodes, i.e., the frequency of the cluster 
motif in those models. Furthermore, the bulky D cell tends to force its surroundings, so 
increases in its density tend to decrease the entropy density and increase the volume added 
per dead-surface (see Sec. Ill A below). Thus, the chemical potential conjugate to ( is the 
most important one for approaching the correct ensemble. 

We can write a more general form for the free energy /(E; () which has an absolute 
maximum at E = 0, ( = Co- (Possibly it might depend on additional density parameters 
also, but for the moment we will stick to this simplest case.) Then the expansion of / around 
this maximum will contain terms such as §i£f(C — Co) 2 ; where is a new elastic constant. 
(Cross terms between E and ( are forbidden by symmetry.) In the present calculations, this 
extra elastic constant is not physically relevant, for reasons which we discuss in Section III D 
below, so we define the entropy as the maximum 5(E) = max^[S(E, ()], making it a function 
only of the phason strain. 



III. TRANSFER MATRIX APPROACH 

In the majority of previous calculations on three-dimensional random tilings, workers 
have made use of Monte Carlo methods to study the elastic properties. In two dimensions 
on the other hand, the transfer matrix method has probably been as important as Monte 
Carlo simulation |7|,|39|-^1[1 . General views of the transfer matrix method as applied to two- 
dimensional quasicrystals may be found in Ref. pifl , Sec. 8.5.1, and in Ref. |J3). Transfer 



matrices provide a potentially exact way to calculate phason elasticity, and are more conve- 
nient for calculating entropies, which can be extracted from a Monte Carlo simulation only 
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after integrating an entropy differential from zero to infinite temperature [fll ,|I4 , 4^] . In the 



case of the twelve-fold 'square-triangle' tiling jO]], a transfer-matrix formulation has made 
possible an exact solution with analytic expressions for the entropy density and the elastic 



constants p4|.|43 



In the case of the canonical cell tiling, the implementation of a Monte Carlo simulation 
has, as noted above, been blocked by the lack of an update move. So instead, we have 
applied the transfer matrix method, using a new technique, which we now describe. 

As described in Section II C, we consider an ensemble of all configurations filling a 'tower' 
with a given, finite base and having periodic boundary conditions in the transverse (x and y) 
directions, but unbounded in the z direction. This reduces our three-dimensional problem 
to one which is essentially one-dimensional, and therefore amenable to a transfer matrix 



treatment. We specialize to the 'maximally random' ensemble of tilings |21| , meaning that 
each tiling configuration has equal weight (the Hamiltonian 7i is zero). 

In Paper I we demonstrated how our towers can be decomposed into layers, and how 
the set of all possible towers may be codified in a 'stacking graph', which represents all the 
different ways of stacking these layers one on top of another. In this section we review this 
development and then discuss how the stacking graph may be turned into a transfer matrix 
whose eigenvalues are related to the free energy of the ensemble of tilings. We will also have 
to deal with a few technicalities associated with the facts that (i) the layers do not have 
equal volume and (ii) we must be able to compute the free energy for any accessible phason 
strain. 



A. Dead-surfaces and the stacking graph 

The basis of our transfer-matrix approach is a rule for representing any given tower of cells 
by a sequence of layers of cells stacked in the z direction, with a one-to-one correspondence 
between possible tilings and sequences of layers. Given that we can formulate a rule for 
dividing our tower of cells into these layers (possible strategies for doing this are discussed 
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below) we can then construct a 'stacking graph' for the complete ensemble of tilings. The 
stacking graph consists of vertices joined by arrows (see Fig. 4 of Paper I, for example) where 
each vertex represents one kind of layer, and an arrow from one vertex to another represents 
a possible way in which the first kind of layer can be followed by the second. Any possible 
structure can then be described by a particular path through the stacking graph. 

There are two different approaches to dividing a system into layers. In the first approach, 
which has been used in all previous transfer matrix studies of random tilings, one slices the 
structure up into slabs of cells, which span the entire cross-section of the tower from side 
to side. The average number of cells in a slab thus scales with the cross-sectional area. 
In two dimensions, the layers can often be mapped onto rows of sites on some lattice. By 
interpreting the stacking direction as a 'time' axis, conditions on the allowed transitions from 
one layer to the next can then be interpreted as local rules for the hopping of a conserved 
set of particles moving in one dimension. Making use of this idea for the case of the square- 
triangle tiling, a Bethe ansatz solution has been found yielding exact analytic values 
for the entropy and elastic constants H3I . 



The decomposition of three-dimensional tilings into slabs is much less promising since 
each layer is itself a 2D random tiling spanning the cross-section of the tower. At every step 
in the construction of the stacking graph, we would need to completely enumerate all the 
possibilities for the next layer, and the number of such possibilities would grow exponentially 
as a function of the cross-sectional area, making the stacking graph extremely convoluted. 
So in Paper I we introduced a different approach, in which towers of canonical cells are 
divided into layers separated from one another by 'dead surfaces' [jHfl . Using this method, 



we generated the stacking graph and applied it to the production of an exhaustive list of 
periodic structures on the square '|' base |24| . 



The notion of a 'dead surface' is linked with that of 'forcing'. Imagine that we build a 
tower of cells, starting from some initial base, by adding nodes one at a time to the open top 
surface. Often, the top surface will have certain crevices which force the addition of a node 
at a particular place. The reason for this is that there are a number of linkages leaving each 

15 



vertex in the surface, and there are only a finite number of different ways in which linkages 
are allowed to meet at a vertex. It is possible therefore, and indeed quite common, for the 
linkages which are already in place to restrict the number of possible choices for the ways 
in which a vertex can be completed. Thus the addition of a linkage or linkages around that 
vertex may be 'forced', and with these forced linkages come new vertices that sit at their 
other ends. These are our forced nodes. 

A surface which has no forced nodes around any of its vertices is called a 'dead surface'. 
We can produce a dead surface by taking a tower of canonical cells and adding all forced 
nodes to the top until no more additions are forced. The only choices we have to make in 
the construction of the tower are the ones we make about which node to add next when 
we are at a dead surface. Thus, a list of the successive dead surfaces and the choices we 
made when we got to them are sufficient to specify the tiling uniquely. The typical number 
of cells between two successive dead surfaces is relatively small, and approaches a small 
limit as the cross-sectional area of the tower diverges, so that the work involved in finding a 
surface does not increase indefinitely with the size of the system. And, most important, the 
stacking graph is sparse, in the sense that there only a small number (usually two) of possible 
successors to a given dead surface (by contrast with the slabs approach, in which there are 
exponentially many as the area of the slabs becomes large) j|5|]. The 'dead surface' and 
'slab' approaches are contrasted in Fig. [3], using the two-dimensional square-triangle tiling 
for pedagogical purposes. 

B. Construction of the transfer matrix 

If we were simply to generate tilings at random by taking random paths through the 
stacking graph and reconstructing the towers of cells to which they correspond, we would 
not generate tilings with the same weights with which they appear in our ensemble. To 
see this, consider a stacking graph in which, say, layer 1 can lead to two different layers, 
numbered 2 or 3; taking a random sequence means we would go to layer 2 or 3 with equal 
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probability. But in the random-tiling ensemble, it may be that layer 2 is more likely than 
layer 3. That will be the case if layer 2 ultimately leads to a larger number of tilings 
afterwards than layer 3 does. If we specify that all tilings should appear with equal weights 
(Ti — 0), that is not the same as saying that all transitions from one layer to another appear 
with equal weights. It is in order to tackle this problem that we introduce the transfer 
matrix. 

As explained in the preceding section, a tiling of M layers is uniquely represented by 
a sequence of dead surfaces and the choices made to continue growing the tiling at each 
surface. We will denote this (a m ,l m ) with m = 1, . . . , M where a m is the label of the 
m th dead surface and l m is the label of the transition a m — > cr m+1 in the stacking graph. 
(Normally there will only be one choice we can make at surface a m that will lead us to 
surface (7 m+1 . But occasionally there will be more than one way to make the transition, and 
in these cases we use the label l m to distinguish between them.) 

We also partition the Hamiltonian into terms which each involve just two surfaces and 
the tiles contained between them: 

M 

H = U(a m , <7 m +i, l m ). (18) 
m=l 

Then the partition function can be written as 

Z= ]T e' n = TiT M , (19) 

{ & m ,1m. } 

where T is the transfer matrix. The general form for the elements of the transfer matrix is 

T aT = J2e- u ^ l \ (20) 
i 

Note that T ar is zero if there is no connection a — > r in the stacking graph; for the particular 
case Ti. — 0, T aT is just the number of ways to get from a to r. 

The total weight in the partition function of all configurations of M layers starting with 
layer a and ending with layer r is then (T M ) CTT . It can be seen that, whatever boundary con- 
ditions we choose for the first and last layers, the free energy per layer in the thermodynamic 
limit is 
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Flayer In ^-max (21) 

where A max is the largest eigenvalue of T. 

We are interested in computing the free energy for a range of phason strains, and also 
for a range of values of the parameter ( (see Section II D). This means we must consider 
(i) how to measure the mean phason strain and density of a given ensemble, and (ii) how 
to generate ensembles with differing phason strain and density. These are dealt with in the 
next two subsections. 



C. Expectations 

Having found the dominant eigenstate of the transfer matrix, we can extract the expec- 
tation of any operator Q that can be written in the form of a sum over the transitions from 
one dead surface to the next: 

M 

, lm) ■ (22) 

m=l 

If v ^ and are the left and right eigenvectors corresponding to A max , then the expectation 
of Q is 

(Q) = E4 i) e- l/(CT ' T '°g(a,r,/)^). (23) 

<T,T,Z 

Thus, for example, the expectation of the Hamiltonian, which is just the average energy per 
layer of the random tiling is 

(U) = vi L h~ u ^U(a, r, l)v[ R \ (24) 

cr,r,/ 

A number of other useful operators besides the Hamiltonian may be written in the 



form (|22|). In particular, having designated a representative node on each dead surface, we 
can define vectors R(<7, r, I) and H(cx, r, I) which are respectively the step in physical space 
and in perp space from one dead surface to the next. Then the total offset in physical and 
perp space are given by 



M 




(25) 



and 



H 



■tot 




(26) 



m=l 



In the thermodynamic limit, the phason strain tensor E is then given by 



H, 



■tot — 



E R, 



■tot- 



(27) 



Six of the nine components of E take fixed values which we already know (see Section II C). 
The remaining three we calculate from (p7|). 

We can also write the total number of nodes as 

M 

Ntot = E N ( a m, Vm+l, l m ). (28) 
m=l 

where N(a m , cr m +i, lm) is the number of nodes added between surfaces m and m+1. Since 
the total volume is given by 



(where a is the length of the edge of the square base), the number density of nodes is given 
by the ratio N tot /V t ot- 



We wish to investigate ensembles of towers of canonical cells with different mean phason 
strains, so as to be able to maximize the random tiling entropy as a function of phason 
strain. In order to vary the average values of the three free phason strain components in 
our ensemble, we need to introduce terms into our Hamiltonian which couple to the phason 
strain, as mentioned briefly in Section II D. For example, if we want a particularly large E xz 
component of phason strain, we should favor those transitions from one surface to another 
which contribute a large shift in the x direction in perp space by comparison with the z 



Vtot — fR-tot]^' 



(29) 



D. Chemical potentials 
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component of the accompanying real space shift. The appropriate form for the terms in 
the Hamiltonian to achieve this is \x ■ A, where A is the vector introduced in Section II C 
composed of the three free strain components, and // is a vector whose components are 
chemical potentials coupling to the phason strain. Then the entropy per layer S\ a y er is given 
by Legendre transformation: 



layer 



In A max = -Slayer + }J> • A (30) 



where A is the mean phason strain. Physically, we are only interested in the case 7i = 0; 
the chemical potentials are added solely as auxiliary fields in order to probe the variation of 
the entropy with phason strain. 

We can similarly introduce a chemical potential /i^ coupling to the the parameter (, or 
equivalently to the packing fraction of nodes. However, we assume that ( is freely varying, 
since we do not envisage placing any physical constraint on our ensemble that would fix 
its value. So the appropriate course of action will always be to choose the value of ^ that 
maximizes the entropy density . 



For a system which did have a non-zero Hamiltonian, the equilibrium value of the internal 
energy (where T = 1) would be given by 

U = S-h-A-kC, (31) 

and the equilibrium free energy, which is also the free energy at which the entropy is a 
maximum, is given by the Legendre transform: 

F = U-S + fi-A + fi c ( = 0. (32) 

So we can either look for the maximum of the entropy, or equivalently we can look for 
the zero of the free energy, to find the equilibrium values of A and (. In our particular 
calculations we maximized the entropy. 

It is not necessary, as it is in some models, that the maximum entropy occur at the 
point at which all the chemical potentials are zero. The zero of the chemical potentials in 
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this model corresponds to the maximum of the entropy per layer, which has no particular 
physical significance. The physically interesting quantity in this case is the entropy per 
unit volume. This quantity is not trivially related to the entropy per layer, since different 
layers make different contributions to the volume of the system. The calculation actually 
performed involves finding the entropy per layer and then converting it to an entropy density 
by dividing by the mean volume per layer within the ensemble: 

S V = §£, (33, 

where R z is the z component of the ensemble average of the operator R(cr, r, I). The mean 
volume per layer will vary with the chemical potentials as we weight layers of different 
volumes more or less heavily. Because of the reflection symmetry of the ensemble in the x 
and y axes, we do in fact expect the minimum entropy density to fall at \i x — // y = 0, but in 
general it will not also fall at /i z = //£ = and so we must probe a range of values of these 
chemical potentials to find it, calculating the entropy density as shown above. 



IV. RESULTS AND DISCUSSION 

Our calculation proceeds as follows. We define the transfer matrix T as in Equation (|19D 
using the list of possible transitions from one dead surface to another for a particular system 
generated from the stacking graph of Paper I. The systems we have looked at are the square- 
base canonical cell tilings with base of length rb, r 2 b, and r 3 b (known as the and '|' 
sizes). Although the transfer matrix quickly becomes large as the system size increases (the 
| system has a transfer matrix some 7000 elements on a side) it is very sparse, so that matrix 
multiplication can be performed quickly, in time of order the rank of the matrix rather than 
the number of elements. The great advantage of the method we have adopted is that the 
quantities we are interested can be calculated by evaluating only the largest eigenvalue of 
the matrix, and its associated left and right eigenvectors. We can find these by simply 
multiplying the matrix many times into a trial eigenvector. Assuming this trial eigenvector 
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has a non-zero component in the direction of the lowest eigenstate of the matrix, this will 
quickly give us the lowest eigenstate, and either one of the left or right eigenvectors of the 
system. We then repeat the process to find the other eigenvector, and from these we can 
evaluate the free energy per layer of our ensemble from Equation fl2~T|) and the free phason 
strain components A from Equation (^7|) as well as the parameter (. Using these we can 
evaluate the entropy density per node from Equation (|33|). 

We then vary the values of the chemical potentials coupling to the phason strain and ( 
to find the maximum of the entropy for the ensemble, and the quadratic variation about 
that maximum to extract the coefficients C\, C2, and C3, and hence the phason elasticities. 
The most straightforward way to accomplish this turns out to be to repeat the calculation 
on a grid of values in the space of the three chemical potentials that make up the vector /x, 
maximizing always with respect to the remaining potential fi^, as discussed in Section III D. 
Then we perform a cubic least-squares fit to the resulting function. The quadratic terms in 
this fit give us our elasticities, and the cubic terms give us an estimate of the error involved 
in calculating second derivatives at the maximum using what is essentially a finite difference 
method. 

For the '|' system, which has a base of length rb and a background phason strain of r~ 5 , 
we find that one of the horizontal components of the phason strain — the one that we call A2— 
can only take one value no matter what sequence of layers we stack one on top of another. 
Thus the corresponding elasticity parameter, C2, is infinite. The other two parameters are 
finite. For the '|' system, which has a base of length r 2 b and a background phason strain of 
r~ 7 , all three C parameters have finite values, but the system is still sufficiently constrained 
that the packing fraction of nodes, or equivalently the parameter (, takes only one value 
independent of the order in which we stack surfaces to construct our tower, so that there 
is no possibility of allowing this parameter to fluctuate and maximizing the entropy with 
respect to it. 

The first system large enough to exhibit the generic elastic behavior typical of towers 
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of canonical cells with large bases is the '|' system, which is also the largest system we 
have tackled. In this system, which has a square base of length r 3 b on each side and a 
background phason strain of r~ 9 , all three components of the strain vector A are free to 
fluctuate under the influence of their corresponding chemical potentials, as is the parameter 
( which measures the density of nodes of the tiling. Figure |2| shows a contour plot of 
the entropy surface max^S^A, ()] as a function of the phason strains Ai and A 3 . The 
maximum entropy density is 6.00 x 10~ 3 , and falls at A\ = and at a small non-zero value 
of A3 (= 0.0077. . .), as predicted. (It also falls at A2 = 0, though this is not evident from 
the figure.) We have also studied the dependence of the entropy on ( and we find that for 
a given E, it depends very sharply on There is a narrow valley in the plot of f(E;Q 
centered around the plane ( = Co(E). In consequence, if we fix (, the apparent dependence 
of /(E) on E is very sharp. (Note that, for the 2/1 and 3/2 cases, ( is in fact a function 
of E — the values of the background phason strain are so high in those cases that the valley 
has zero width.) 

If we just omit (set to zero) the fi^ parameter, thereby accepting the ensemble which 
maximizes the entropy per layer at a given phason strain, then we find an entropy density 
which differs from the actual value by only about three parts in 10 5 . We assume that the 
difference produced by including terms in the Hamiltonian coupling to other densities of 
local tiling patterns is smaller still, and therefore that we were justified in neglecting them 
in this calculation. 

The value of the parameter £ at the maximum of the entropy for the f ensemble was 



0.334, which is close to the "magic" value of ( = 0.317 obtained at the end of Ref. [25 



The corresponding packing fraction <fi for spheres of radius c placed on the nodes may be 



determined from formulas given in Sec. Ill A of Ref. [ 2l| . After converting from the length 
units used there, we obtain 

, vrv/3 



16 
where 
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7^ + 2(1 -, f )(i-§; 



(34) 



7] = -(1 -r 9 detE) (35) 

and detE = I? 2 ^ -* in our case. For the | ensemble, we find a value of <p — 0.606. Values 
for the other system sizes are given in the table. 

The packing fraction at which the maximum entropy occurs can be written roughly as 

0(E) = 0.6064 - 0.026(A 3 -4 0) ) (36) 

for the | case, with virtually no dependence on A\ or A^. (It seems plausible that, with 
more general values of phason strain, it would depend on the determinant det E of the strain 
tensor.) In other words, the optimal packing fraction is practically constant (varying ~ 1% 
over the range of possible phason strains.) In the \ and | cases, is a unique function of 
the phason strain. In these systems the coefficient of variation with A 3 is respectively 3.5 
and 0.5 times as large as in Eq. (|3"6"D. 

Values for the elastic parameters C\, C2, C3, and for the phason strain at equilibrium 
Af } for the various systems are given in Table I. The values do not fit the relations fllBD 
very well, presumably because the system size is small. However, it is possible to make out 
some trends in the data as the system size increases: 

(i) There is a clear tendency to oscillate between high and low values of the elastic con- 
stants from one approximant to the next. 

(ii) superimposed on this oscillation, there is also a strong tendency for C m to decrease 
with increasing system size. 

The first of these trends is to be expected if the errors are caused by the large background 
phason strain B, in view of Equation ([13]). The same oscillation is seen in the Monte Carlo 



simulations of the rhombohedral tiling by Shaw et al. |I3" |. (See, for instance, Figure 2 in 
that paper.) 

We suggest that (ii) is also a finite- B effect. The phason strain tensor components define 
a 9-dimensional space, and there is a domain in this space (centered on zero) of allowed 
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phason strains. The entropy density vanishes with singular derivatives at the boundary of 
this domain. When B is large, the phason strain tensor is necessarily closer to the boundary 
of the allowed domain, and the second derivatives of the entropy density (including the 
coefficients C m ) are consequently larger. In addition, the tendency for Sy to be larger 
for smaller systems also enhances C m , because increasing the maximum value of a function, 
while at the same time decreasing the interval over which it must rise from to its maximum 
and fall back to zero, obviously forces it to have a larger (negative) second derivative. 

So what can we say about the extrapolation of the data to infinite system sizes? The 
entropy density Sy at least, appears to be far from its asymptotic value in the systems we 
have studied here. Indeed, on the basis of these data alone, one could not really rule out 



Sy = in the infinite-system limit [47|] ! 



For the coefficients C m , in view of the large values of C 2 for the 2/1 and C3 for the 3/2, 
we have extrapolated the inverses C" 1 . We have assumed the behavior 

1 1 



+ XB, (37) 



where Cffl is the value of C m measured for the n th approximant and B is background phason 
strain. For each C m , we have extrapolated from the pair (2/1,5/3) and the pair (3/2,5/3) 
and taken the mean of these as the extrapolated value, writing error bars to span the two 
independent extrapolations ||48|| . The results are given in the "extrap" column of the Table. 

We have also performed a slightly different analysis of the 5/3 system using just the 316 
approximants described in Paper I. This calculation is detailed in Appendix A, and in the 
column marked "| (F) in the table. 

Each of these different estimates give us four quantities C\, C 2 , C3, and A3 , but there 
are only two independent parameters Ki, K 2 to fit to them. This places constraints on the 
four quantities with which their measured values are not entirely consistent. However, we 
can certainly conclude that K\ « 1, and that K 2 is negative and fairly large (possibly close 
to the instability limit ([TDD ) . This sign of K 2 is opposite to that in the rhombohedral tiling, 
a fact that we can explain qualitatively by the following argument. 
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From ([16]) we have 



K, = ^| (38) 



so that K2 < if and only if ft, x fluctuates more than h y as one walks along a path in the z 
direction. Examining the tiling, we find that the linkages along such a path are dominated 
by the 2-fold b linkage that runs entirely along z, and by the 3-fold c linkage (±r, 0, r 3 ). In 
our basis (§J), these c linkages have perp space components only in the x directions, making 
the x-fluctuations large compared with the y ones, which in turn makes K2 negative. On the 
other hand, the linkage vectors in the rhombohedral tiling are just the basis vectors in @; 
thus the dominant steps in the direction of the path are 77(0, ±1, r) which has fluctuations 



in the h y direction rather than the h x direction. Thus for that tiling K 2 > [|T2| , p~3 . 

For applications of random tilings to structural fitting, the parameter of greatest interest 
is the effective phason Debye- Waller factor B -1 , which corresponds to the extra variance in 
perp space due to the randomness of the tiling. The experimental value of B -1 is of the order 
of IrA 2 . In Appendix B, using the ensemble of 316 packings we find a crude theoretical 
estimate of B -1 for the canonical-cell random tiling which is somewhat smaller. 

Finally, we would like to consider whether there is the possibility of performing a cal- 
culation for a larger system. Probably there is. It would tax the power of the available 
computing resources, but it should be possible. We estimate (see Paper I) that the next 
largest size of system — the | system — should have a stacking graph of about one million 
dead surfaces. Given that the average number of nodes added between one dead surface 
and another tends to a constant as the system becomes large, we expect the time taken to 
construct the graph to scale as the number of dead surfaces, and so we would require about 
200 times as much CPU time to perform the calculation as we did for the | system. With 
today's high-performance computing resources, such a calculation would be just within the 
bounds of possibility. 
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V. CONCLUSIONS 



Using the results of a previously-presented method for breaking down and cataloging 
three-dimensional random tilings, we have defined a transfer matrix whose eigenstates de- 
scribe the properties of an ensemble of random-tiling configurations in the shape of towers 
with a finite area base and infinite extent in the z direction. The ensemble has a mean 
phason-strain and packing fraction controlled by chemical potentials whose values we can 
choose. The transfer matrix is sparse and so can be efficiently multiplied into a trial eigen- 
vector to find the largest eigenvalue. This gives us the free energy density and thus the 
entropy density for the ensemble. Maximizing this with respect to the components of pha- 
son strain and the packing fraction we can find the equilibrium values of these quantities, 
and the variation of the entropy density about the maximum, which defines the phason elas- 
ticity constants which are related to experimentally measurable quantities such as phason 
Debye- Waller factors, and diffuse scattering close to Bragg peaks. 

We have applied our method to the 'canonical-cell tiling' p5[ for towers with a square 
base with sides taking the three smallest non-trivial lengths possible. For the two smaller 
of these systems, we find that the systems have additional constraints making their elastic 
behavior different from the generic behavior expected of large systems. The largest system 
we have studied, which we refer to as the '|' system is the smallest one that exhibits the 
same elastic behavior as a large system. From this system we have extracted a value of 
6.00 x 10 -3 for the entropy per unit volume. The two fundamental elastic constants, K\ and 
K%, are over-determined by the transfer matrix calculation because elasticities corresponding 
to the strains in various directions are related by symmetries of the system. This gives us 
a way to estimate the (finite-size) errors in our calculation, and it turns out that the values 
for K% and K 2 are not very accurate for this smallest of systems. Our best estimate is that 
Ki lies in the range 0.8 to 1.1, and K 2 < 0. In addition, the requirement that the tiling be 
elastically stable with respect to other phases means that K2/K1 > — |. 

It is interesting to compare our results for the canonical-cell tiling with those for the 
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rhombohedral random tiling [pT2] , p~3|| . In doing so, we must be careful to allow for the different 
lengths of the linkages making up these tilings. 

In the rhombohedral tiling [14]] the entropy per node Sn is ~ 0.24 ± 0.02, whereas in 



the canonical-cell tiling it is only ~ 0.047 ± 0.01. This is a result of the fact that the 
canonical-cell tiling is more constrained in the states it can take than the rhombohedral tiling. 
This constraint is also responsible for limiting the canonical-cell tiling to only 32 distinct 
local node environments ipER . By contrast, the node environments of the rhombohedral 
random tiling are chosen from 10 527 of different possibilities ||50|| . 

To compare the elastic constants of the two tilings, we first define a characteristic unit of 
phason strain E to be the ratio of the perp- and real-space displacements of the fundamental 
linkage. E is 1 and r -3 for the rhombohedral and canonical-cell tilings, respectively. Then 
a dimensionless measure of the elastic constant K\ is the corresponding decrease in free 
energy per node for a phason strain of magnitude Eq\ 

AF N = X -K x Eljn, (39) 

where n is the packing fraction of nodes. For the rhombohedral tiling, n = 0.6498, K\ = 
0.81 ± 0.01 and K 2 = 0.495 ± 0.02, hence AF N = 0.263; for the canonical-cell tiling, using 



the extrapolated values AFn = 0.28, very similar to the rhombohedral case [[51 . 

Although the results for the elastic properties of the small canonical-cell systems show 
signs of being strongly influenced by finite-size effects, we believe that it should be feasible 
using modern supercomputing resources to perform a similar calculation for a larger system 
which should yield a more accurate estimate of the elastic constants for the infinite random 
canonical-cell tiling. We also believe that the method we have presented, which is fundamen- 
tally different from transfer matrix methods previously used in the study of random tilings, 
should be applicable to any random-tiling model in any number of dimensions. Even if the 
canonical-cell model does not ultimately turn out to be as good a model of real icosahedral 
phases as some other (yet to be proposed) tiling, we believe this method will prove useful in 
the calculation of experimentally measurable properties of quasicrystals. 
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APPENDIX A: FOURIER MODE ANALYSIS 

In the past, elastic constants have been extracted from Monte Carlo simulations of finite 
random-tilings with periodic boundary conditions, through measurement of the equilibrium 
fluctuations of the phason strain. We have carried out such an analysis for the present 
system, to provide independent (albeit inferior) estimates of the elastic constants. In lieu 
of a simulation, we have used the exact expectations for the ensemble with 316 microstates 
computed in Paper I, which are packings of a cube with side L = r 3 b. 

The elastic free energy (|8|) can be expanded around a state of zero phason strain and 
rewritten as a sum over Fourier modes, 

F tot = 5>(-q)K(q)h(q) (Al) 
q 

where h(q) is the Fourier transform of the height field (gj) normalized as in Ref. | 13| . In the 
basis of (|]) and (f|), the stiffness coefficients are |30| 

ffy(q) = [KM? - ^ 2 [(i|q| 2 + r~ l q 2 i+l - rg 2 +2 )]^- - 2^(1 - <%) . (A2) 

It then follows that 

(Mq)Mq)> = {k^W. (A3) 

In order for h(r) to be well defined in a cell with periodic boundary conditions, we must 
adopt rationally- related perp space vectors as used in |I5]| . In the present case, the perp 
space basis vector i](r, 1, 0) is replaced by ?/(5, 3, 0) where 7/ = (5 2 + 3 2 ) -1 / 2 , and similarly 
for the other vectors in Equation (^J). 
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For our Fourier transform we used the crude definition 

v 1 / 2 



h(q) 



N 



(A4) 



following Ref. (Interpolations of h(r) between tile vertices were used by Refs. (15 

and 0.) 

We report results only for the smallest wavevector q = (0, 0, 2tc/L), since even this is too 
large to truly be the long-wavelength limit. For that q value, K only has diagonal elements 
and they are simply Ku = |q| 2 Cj with Cj given by Equation (|l^) |52| |. Finally we obtain 
the estimates Ci = |q| _2 (|h(q)| 2 ) -1 which were reported in the table. (The results from the 
next larger q value, (2tt/L,2tt/L,0), are also consistent with Ki < 0, and K2/K1 being of 
order unity, but with K\ and K 2 each increased by a factor of at least 2.) 



APPENDIX B: THE PHASON DEBYE- WALLER FACTOR 

To define the perp-space Debye- Waller factor, one assumes that there exists an ideal, 
perfectly quasiperiodic structure made of canonical cells and that our random tiling (when 
represented as a surface in 6-dimensional space) differs from this ideal structure by random 
displacements in the perp-space direction, which have a Gaussian distribution with each 
component having variance 2B^. It follows that each actual structure factor is reduced 
from its ideal value by the factor exp(— B^IG -1 ! 2 ), where G -1 is the perp-space component 
of the reciprocal lattice vector. 

Fitting of experimental data found B ± = 0.39rA 2 for i(AlCuLi) |3f and B ± = 0.70rA 2 
for i(AlPdMn) |54] (in units where the basis vectors in Eqs. ([I]) and (0) have length = 
5.1rA and 4.65rA respectively.) In the dimensionless units used in this paper, w 0.015 
and 0.034, respectively. If the variance of the microscopic coordinates (0) of the nodes in 
the ideal structure is (|r J -| 2 ) , and in the actual structure their variance is ( | r — 1 2 ) , then 

(|r ± | 2 )-(|r ± | 2 ) = 6B^ (Bl) 

exactly. 
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Using Eq. (|B1|), we have estimated B -1 from the ensemble of 316 packings of the cube 
of side r 3 b (the so-called '|' approximants) constructed in Paper I. Since the quasiperiodic 
ideal structure is unknown, (|r | 2 } is unknown. However, the packing in our ensemble 
with minimum perp-space variance ((r^l 2 ),^ has a high symmetry and is probably a true 
approximant of the quasiperiodic structure; the perp-space variance of this structure should 
not be too different from that of the ideal infinite structure. We find (|r J -| 2 ) miri = 0.2154 
and (Ir^l 2 ) = 0.2758. Hence 

« ^(Ir^l 2 ) - (\r x \ 2 ) min « 0.010 (B2) 
o 

This is significantly smaller than the experimental value reported for i(AlCuLi), and much 
smaller than that for i(AlPdMn). 
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TABLES 



size 


edge 


B 


Ci 


c 2 


c 3 


AO) 


<p 


Sy 


2/1 (TM) 


rb 


0.090 


0.890 


OO 


1.282 


0.0185 


0.6004 


0.0292 


3/2 (TM) 


T 2 b 


-0.034 


2.101 


5.835 


10.14 


0590 


0.6051 


0.0103 


5/3 (TM) 


T% 


0.013 


0.452 


2.302 


0.494 


0.0077 


0.6064 


0.0060 


extrap. 


OO 





0.48 ± 0.08 


1.54 ±0.35 


0.48 ±0.10 


0.014 
±0.008 


0.607 
±0.001 


0.0041 
±0.0007 


5/3 (F) 


T 3 b 


0.013 


0.205 


2.663 


0.645 


0.013 


0.606 


0.0036 


formula 


T n~2 b 


r -(2n+3) 


K x + 1.285K 2 


Ki - 0.951K 2 


Ki - 0.333K 2 


-2BK 2 






Ki -0.333^2 







TABLE I. Results for the elasticity coefficients and equilibrium phason strain of various sizes 
of system. The rows labeled TM were calculated from the transfer matrix method, and the row 
marked MC was calculated by analyzing the 316 cubic | approximants from Paper I using the 
method employed by Shaw et al. in their Monte Carlo studies of the square-triangle tiling (see 
Appendix A). The "extrap." row is an extrapolation of the transfer matrix results to a system of 
infinite size. 
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FIGURES 




FIG. 1. Dividing a strip of square-triangle tiling into layers, (a) 'Slab' approach (as used 
by 0|). (b) 'Dead-surfaces' approach (as used in the present paper for the 3D canonical-cell 
tiling. 
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FIG. 2. Contour plot of the entropy as a function of the phason strain components A± and A% 
in the region of the maximum entropy. 
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